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Detailed protocol 


We modeled the long-term dynamics and diversity of ecological communities using the well-known generalized Lotka-Volterra (gLV) model, | 


modified to include dispersal from a species pool: 
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where Nj is the abundance of species / (normalized to its carrying capacity), aj is the interaction strength that captures how strongly species / inhibits the 
growth of species / (with self-regulation aj; = 1), and D is the dispersal rate from an outside species pool to the focal community. For simplicity and without 


qualitatively changing our results, we considered the same growth rate 
7 
= 1 and the same carrying capacity 


Ki 


= 11 for all species in the main text. Fig. S6 shows that sampling growth rates from a uniform distribution has little effect on the phase diagram of survival 
fraction and fluctuation fraction. Fig. S6 shows that sampling carrying capacities from a normal distribution increases the partial coexistence phase while 
shrinking both the full coexistence phase and fluctuation phase, but does not affect the order of the phases. We further tested the theoretical predictions when 
considering the existence of positive (facilitative) interspecies interactions (Fig. S9) and varying the symmetry of the interaction matrix (Fig. S8). We also 
considered different dispersal rates (Fig. S7), and the effects of incorporating daily dilutions (Fig. S9) in thesein silico communities. These additional results 
show that our qualitative phase diagrams and conclusions are robust to different choices of ecological network structure and parameters. Although the 
patterns of ecological diversity and dynamics do not change as the dispersal rate varies from D=107 to D=10% (Fig. S7), we found that communities with 
zero dispersal rate exhibit lower fluctuation fraction and survival faction in the persistent fluctuation phase. Our results show that non-zero dispersal rates can 
sustain persistent fluctuations. 


To show that symmetry and anti-symmetry do not qualitatively change the phase diagram, we simulated communities in two scenarios of non-zero 
reciprocity, 


y=0g 
and 
y=-05 
, where the reciprocity of interactions is given by 
. The qualitative patterns and order of transitions in the phase diagram (Fig. S8) are robust to the presence of reciprocity, although 
¥ 
shifts the stability boundary (solid line in Fig. S8). Positive and negative reciprocity decreases and increases the values of S and < 
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> at which communities lose stability, respectively. Moreover, positive (negative) reciprocity yields lower (higher) survival fraction and fluctuation fraction of 
communities. At full symmetry and anti-symmetry ( 
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=1,-1) there is no fluctuating phase (27), but those do not appear to be relevant from the experiment pair-competition results (Table. S1-S3). 


To test whether the existence of positive (facilitative) interactions in the ecological network will change our conclusions, we sampled values of 
kul 
from a uniform distribution [- 
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bed 
], where 

kad 
varies between [0, 1.4] on the phase diagram. In this simulation, the linear interaction function in the gLV ( 
) is replaced with Monod function ( 


) to avoid unbounded growth due to positive interactions (27, 42). We observed similar patterns of survival fraction and fluctuation fraction between pure 
competitive interactions and considering positive interactions (Fig. S9). The first and second moment of the distribution of 


be] 

should be considered to quantify interaction strength in the stability criteria (17). Here we use Std ( 
be] 

) to quantify the interaction strength because the first moment of 


bul 


distribution is zero. Our results demonstrate that the existence of three phases (full coexistence, partial coexistence, persistence fluctuation) and the order of 
transitions are robust to the interaction types in the model. 


All simulations used the Runge-Kutta method on Matlab to numerically solve the LV equations (with an integration step of 0.05). A definition of 100 


i 


100 pixels was used for each phase diagram, linearly segmenting the parameter space in the ranges <aj> © [0, 1.5] and S © [1, 100]. In each phase 


diagram, each pixel shows the average result for 103 simulations. The total simulation time is 104 to guarantee the survival fraction and fluctuation fraction 
have reached steady states as shown in Fig. S5. 


To test whether the total biomass fluctuation is consistent with species abundance fluctuation in our simulations, we simulated the time series of community 
biomass under various conditions (Fig. S3). To quantify the dynamics of biomass in the simulation, we calculated the sum of species abundance ( 


). Our results demonstrate the fluctuation in species abundance is in agreement with fluctuation in total biomass. 

The similar fluctuations between replicates in some experiments (e.g., Fig. 2, medium nutrients, S=48) could be explained by the slow divergence that 
chaotic trajectories can exhibit during moderately long-time windows, or, alternatively, by possible limit cycle oscillations (Fig. S3). We focused on chaotic 
fluctuations when discussing the model predictions because previous theory shows that all persistent fluctuations will be chaotic as number of species in the 
pool S grows (27), though limit cycle oscillations only happen under finite S. 


Below are the Matlab code: 


clear;clc;close all; 

r_mean=1.0; % mean growth rate 

S=1:1:100; % number of species 

DD=1e-6; % migration rate 

A_mean_range=[0.01 1.5]; % mean inhibition alpha 
ystep=100; %for A_mean 

T_real=10000; %total evolution time 

step=0.05; % time interval 

T=ceil(T_real/step); Yonumber of circulation 
num_sim=10000; % number of simulation 

period=2000; % duration for calculating richness 

% composition_full=zeros(length(S)*ystep*num_sim,length(S)); % record the final compositions of all simulations 
abundance_threshold=1e-3; %threshold for survival species 
fluc_record=zeros(length(S)*ystep,num_sim); 
diversity_record=zeros(length(S)*ystep,num_sim); 


A_mean=linspace(A_mean_range(1),A_mean_range(2),ystep); 
for mm=1:length(S) 
for pp=1:ystep 
for hh=1:num_sim 


composition= LV_compute(r_mean, S(mm), DD, A_mean(pp),T,step); 
diversity_index=max(composition(T-period+1:T,:)); richness=find(diversity_index>abundance_threshold); 
diversity_record((mm-1)*ystep+pp,hh)=length(richness)/S(mm); 
fluc_CV=std(composition(T-1e5:T,:))./mean(composition(T-1e85:T,:)); fluc_CV(isnan(fluc_CV))=0; 
fluc_record((mm-1)*ystept+pp,hh)=mean(fluc_CV(fluc_CV>0)); 
end 
end 
end 


function [composition] = LV_compute(r_mean, S, D, A, T, step) 


N_mean=1/(S*A*sqrt(pi/2)); 

sigma3=N_mean/12; 

NN=zeros(T,S); 

r=ones(S,1); 

% Interaction matrix 

% AA=exprnd(A,S,S); 
AA=rand(S,S)*A*2; 
AA=AA-diag(diag(AA))+eye(S); 

% considering carrying capacity 

% K=normrnd(1,sqrt(1/12),[S,1]); 

% for i=1:S 

%  forj=1:S 

% AAI j)=AA(iJ)*KG)/K(i); 

% end 

% end 

% Initial composition 
NO=abs(normrnd(N_mean,sigma3,1,S)); 
NN(1,:)=NO; 

% Migration from source 
Migra=ones(1,S)*D; 


% compute the dynamics in each patch 
for i=2:T 


for j=1:S 
k1=r(j)*NN(i-1,j)*(1-AAG,:)*(NN(i-1,:)'))*step; 


k2=r(j)*(NN(i-1 j)+k1/2)*(1-AAG,:)*(NN(i-1 ,:)))*step: 
k3=r(j)*(NN(i-1 j)+k2/2)*(1-AA(:)*(NN(i-1 ,:)))*step: 
k4=r(j)*(NN(-1 j)+k3)*(1-AA(:)*(NN(i-1 ,:)))*step: 


NN(i,j)=NN(i-1,j) +(1/6)*(k1+2*k2+2*k3+k4); 


end 
NN(i,:)=NN(i,:)+Migra; 


end 
composition=NN; 
end 


clear; clc; close all; 

load ./uniform3.mat; 

load feasible_exp.mat 

x_step=length(S); y_step=length(A_mean); 

miu1=A_mean(1); miu2=A_mean(length(A_mean)); S1=S(1); S2=S(length(S)); 
color12=[0.7 0.7 0.7]; 

phase_indicate=mean(diversity_record,2); fluc_digi=fluc_record; 
phasenew_indicate=zeros(length(phase_indicate),1); 


critical_CV=0.001; 
for i=1 :length(phasenew_indicate) 
for j=1:num_sim 
if fluc_record(i,j)<critical_CV 
fluc_digi(i,j)=0; 
else 
fluc_digi(i,j)/=1; 


end 
end 
end 


fluc_indicate=mean(fluc_digi,2); 


phase_std=std(diversity_record,0,2)./phase_indicate; 

phasecolor=reshape(phase_indicate,[y_step,x_step]); phasecolor_std=reshape(phase_std,[y_step,x_step]); 
phasecolor_fluc=reshape(fluc_indicate,[y_step,x_step]); phasecolor_new=reshape(phasenew_indicate,[y_step,x_step]); 
phasecolor_new1=zeros(y_step,x_step); 

x2=linspace(miu1 ,miu2,y_step)'; x1=sqrt(linspace(S1,S2,x_step))'; x2(1)=0; 

X1=repmat(x1'y_step,1); X2=repmat(x2,1,x_step); 


figure(1) 
s=pcolor(X1.42,X2,phasecolor);set(s, 'LineStyle’,'none'); colorbar; hold on; 
caxis([0.0 1.1]) 
plot(xx3, yy3,'k-’,'linewidth’,5,'color’,color12) 
plot(ss,fea,'k-.',‘linewidth’',5,'color’,color12); 
axis([S1 S2 0 1.5]); set(gca,'FontSize',28); 
xlabel('Size of species pool, {\itS}','FontSize',30); ylabel('Interaction strength, <\it\alpha_i_j>','FontSize',30); 
xticks([1 20 40 60 80 100)); yticks({0.0 0.5 1.0 1.5]); box on 


figure(2) 
p=0.9; 
sp=csaps({x1.2,x2},phasecolor_fluc,p); 
v=fnval(sp,{x1.42,x2}); 
s=pcolor(X1.*2,X2,v);set(s, ‘LineStyle','none'); colorbar; hold on; 
caxis([0.0 1]) 
plot(xx3, yy3,'k-','linewidth',5,'color'’,color12) 
plot(ss,fea,'k-.',‘linewidth’,5,'color’,color12); 
axis([S1 S2 0 1.5]); set(gca,'FontSize',28); 
xlabel('Size of species pool, {\itS}','FontSize',30); ylabel('Interaction strength, <\it\alpha_i_j>','FontSize',30); 
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